Let’s import at all of the data that Simon pulled:

library(tidyverse)
library(gridExtra)
library(maps)
library(ggthemes)

DeltasClean <- read_csv("../data/out/deltas_clean_v2.csv") 
Parsed with column specification:
cols(
  Delta = col_character(),
  location = col_character(),
  surface = col_character(),
  year = col_double(),
  month = col_double(),
  ndvi = col_double(),
  red = col_double(),
  evi = col_double(),
  savi = col_double(),
  gr = col_double(),
  ndssi = col_double()
)
DeltaLocations <- read_csv("../data/DeltaLocations.csv")
Parsed with column specification:
cols(
  Deltas = col_character(),
  Lat = col_double(),
  Lon = col_double()
)

As a reminder, for each of the 47 deltas there are measurements of Land & Water areas at Upstream, Downstream and ‘Middle’ locations on the delta. We first lump all the observations together, and look to see which Deltas have many observations:

#counts per delta
DeltaCounts <- count(DeltasClean, Delta)
DeltaCounts

Now, by each month.. where the colorbar represents the number of observations (n) for each month for a given delta:

DeltaObsPerMonth <- count(DeltasClean, Delta, month)

ggplot(DeltaObsPerMonth, aes(y = Delta, x = month, fill=n)) + geom_tile() + 
  scale_x_discrete(limits = c(1:12), breaks = c(1:12)) +
  expand_limits(x = c(1,12)) + 
  scale_fill_gradient( trans = 'log' )

In the above heat map, dark colors (and no color) represent data paucity (and data gaps). Deltas with light colors (e.g., the Parana, Nile, Ebro, Colorado, Brahmani) have lots of data, spread out through the months of the year.

I’ll remove/subset the deltas with sparse coverage (specifically, months with no coverage)….


# need 10 data points per month for NDSSI and NDVI
EnoughObsPerMonth <- DeltasClean %>% ungroup() %>%
  count(Delta, month, surface) %>% 
  group_by(surface) %>%
  filter( n >= 5)

#find deltas missing a given month of observations
DeltaMonthCounts <- EnoughObsPerMonth %>%
  ungroup() %>%
  count(Delta)

# need 12 months of water and land obs, so 24 mo total
EnoughMonths <- DeltaMonthCounts %>%
 filter( n == 24)

CompleteObsDeltas <- pull(EnoughMonths, Delta)

#remove them
DeltasCleaner <- DeltasClean %>%
  filter(Delta %in% CompleteObsDeltas)

#add the real dates in month date format
DeltasCleaner$date <- as.Date(paste(DeltasCleaner$year, DeltasCleaner$month, "01", sep="-"), "%Y-%m-%d")

#remove intermediate data
rm(CompleteObsDeltas, EnoughMonths, DeltaMonthCounts)

#EnoughMonths

and extract some metrics; specifically I will make a timeseries of NDVI and NDSSI for each delta using the mean value for each month.

#take the mean NDVI and NDSSI for each month, for each delta
DeltaMeans <- DeltasCleaner %>%
  group_by(Delta, month, surface) %>%
  summarize(MeanNDVI = mean(ndvi, na.rm = TRUE), MeanNDSSI = mean(ndssi, na.rm = TRUE))

#make a 9 column data frame with:
#delta, 
#max and min NDVI month, 
#NDSSI max and min month, 
#max and min values for both NDVI and NDSSI

#####

DeltaMaxNDVI <- 
  DeltaMeans %>% 
  filter(surface == 'Land')  %>% 
  select (-c(MeanNDSSI)) %>% 
  group_by(Delta) %>% 
  slice(which.max(MeanNDVI)) %>% 
  rename(MaxMeanNDVImonth = month, MaxMeanNDVI = MeanNDVI)

DeltaMaxNDSSI <- 
  DeltaMeans %>% 
  filter(surface == 'Water')  %>% 
  select (-c(MeanNDVI)) %>% 
  group_by(Delta) %>% 
  slice(which.max(MeanNDSSI)) %>% 
  rename(MaxMeanNDSSImonth = month, MaxMeanNDSSI = MeanNDSSI)

DeltaMinNDVI <- 
  DeltaMeans %>% 
  filter(surface == 'Land')  %>% 
  select (-c(MeanNDSSI)) %>% 
  group_by(Delta) %>% 
  slice(which.min(MeanNDVI)) %>% 
  rename(MinMeanNDVImonth = month, MinMeanNDVI = MeanNDVI)

DeltaMinNDSSI <- 
  DeltaMeans %>% 
  filter(surface == 'Water')  %>% 
  select (-c(MeanNDVI)) %>% 
  group_by(Delta) %>% 
  slice(which.min(MeanNDSSI)) %>% 
  rename(MinMeanNDSSImonth = month, MinMeanNDSSI = MeanNDSSI)


#join into 1 dataframe
DeltaMaxMin <- left_join(DeltaMaxNDVI, DeltaMaxNDSSI, by = 'Delta') %>% 
  left_join(.,DeltaMinNDVI, by = 'Delta') %>% 
  left_join(.,DeltaMinNDSSI, by = 'Delta')

#remove intermediate data
rm(DeltaMaxNDVI, DeltaMaxNDSSI, DeltaMinNDSSI,DeltaMinNDVI)

And now we can look at the phase shifts between these two time series (the timeseries of mean NDVI and mean NDSSI). Here are the phase shifts (in month) for each delta:

#compare offset
DeltaMaxMin <- mutate(DeltaMaxMin, 
                      MinOffset = if_else(abs(MinMeanNDVImonth - MinMeanNDSSImonth) > 6, 
                                          12 - abs(MinMeanNDVImonth - MinMeanNDSSImonth),
                                          abs(MinMeanNDVImonth - MinMeanNDSSImonth)),
                      MaxOffset = if_else(abs(MaxMeanNDVImonth - MaxMeanNDSSImonth) > 6, 
                                         12 - abs(MaxMeanNDVImonth - MaxMeanNDSSImonth),
                                          abs(MaxMeanNDVImonth - MaxMeanNDSSImonth)),
                      OffsetDiff = abs(MaxOffset - MinOffset),
                      rangeNDVI = (MaxMeanNDVI - MinMeanNDVI), 
                      rangeNDSSI = (MaxMeanNDSSI - MinMeanNDSSI),
                      halfPeriodNDVI = if_else(abs(MaxMeanNDVImonth - MinMeanNDVImonth) > 6, 
                                          12 - abs(MaxMeanNDVImonth - MinMeanNDVImonth),
                                          abs(MaxMeanNDVImonth - MinMeanNDVImonth)),
                      halfPeriodNDSSI = if_else(abs(MaxMeanNDSSImonth - MinMeanNDSSImonth) > 6, 
                                          12 - abs(MaxMeanNDSSImonth - MinMeanNDSSImonth),
                                          abs(MaxMeanNDSSImonth - MinMeanNDSSImonth)), )

# DeltaMaxMin <- 
#   DeltaMaxMin   %>%
#   select (c(Delta, MinOffset, MaxOffset, OffsetDiff, rangeNDVI, rangeNDSSI,MaxMeanNDSSI,MinMeanNDSSI,MaxMeanNDVI,MinMeanNDVI)) 

DeltaMaxMin

ggplot(DeltaMaxMin, aes(y = Delta, x = MaxOffset)) + geom_point() + 
  scale_x_discrete(limits = c(1:6), breaks = c(1:6)) +
  expand_limits(x = c(0,6))  + 
  ggtitle("MaxOffset")


ggplot(DeltaMaxMin, aes(y = Delta, x = MinOffset)) + geom_point() + 
  scale_x_discrete(limits = c(1:6), breaks = c(1:6)) +
  expand_limits(x = c(0,6))  + 
  ggtitle("MinOffset")


ggplot(DeltaMaxMin, aes(y = Delta, x = OffsetDiff)) + geom_point() + 
  scale_x_discrete(limits = c(1:6), breaks = c(1:6)) +
  expand_limits(x = c(0,6))  + 
  ggtitle("Offset Difference")

Now let’s examine the histograms of all 31 deltas… The months with the greatest mean NDVI, months with gretaest mean NDSSI, the monthly offset, and the skew of the NDSSI and NDVI timeseries.

ggplot(DeltaMaxMin, aes(x = MaxMeanNDVImonth)) + 
  geom_dotplot(binwidth = 1,dotsize = 0.5) + 
  scale_y_continuous(NULL, breaks = NULL) + 
  scale_x_discrete(limits = c(1:12), breaks = c(1:12)) + 
  labs(x = "Month") +
  ggtitle("Month of maximum mean NDVI")


ggplot(DeltaMaxMin, aes(x = MaxMeanNDSSImonth)) + 
  geom_dotplot(binwidth = 1,dotsize = 0.5) + 
  scale_y_continuous(NULL, breaks = NULL) + 
  scale_x_discrete(limits = c(1:12), breaks = c(1:12)) + 
  labs(x = "Month") +
  ggtitle("Month of maximum mean NDSSI")


ggplot(DeltaMaxMin, aes(x = MaxOffset)) + 
  geom_dotplot(binwidth = 1,dotsize = 0.25) + 
  scale_y_continuous(NULL, breaks = NULL) + 
  scale_x_discrete(limits = c(0:6), breaks = c(0:6)) + 
  labs(x = "Months") +
  ggtitle("Months Offset between NDVI and NDSSI")


ggplot(DeltaMaxMin, aes(x = halfPeriodNDVI)) + 
  geom_dotplot(binwidth = 1,dotsize = 0.25) + 
  scale_y_continuous(NULL, breaks = NULL) + 
  scale_x_discrete(limits = c(0:6), breaks = c(0:6)) + 
  labs(x = "Months") +
  ggtitle("half period NDVI ")


ggplot(DeltaMaxMin, aes(x = halfPeriodNDSSI)) + 
  geom_dotplot(binwidth = 1,dotsize = 0.25) + 
  scale_y_continuous(NULL, breaks = NULL) + 
  scale_x_discrete(limits = c(0:6), breaks = c(0:6)) + 
  labs(x = "Months") +
  ggtitle("half period NDSSI")

Ok, so the idea is that peak NDSSI is more effective if it occurs at moderate NDVI, so let’s look at the NDVI value for the months with peak NDSSI.

#extract NDVI value for each delta a the month of max NDSSI value

DeltaNDVIatMaxNDSSI <- DeltaMaxMin %>%
  select(Delta,MaxMeanNDSSImonth)

DeltaMeansToJoin <- DeltaMeans %>%
  filter(surface == 'Land')

DeltaNDVIatMaxNDSSI <- left_join(DeltaNDVIatMaxNDSSI, DeltaMeansToJoin, 
                                 by = c('Delta', 'MaxMeanNDSSImonth' ='month'))
 
DeltaNDVIatMaxNDSSI <- DeltaNDVIatMaxNDSSI %>%
  select (-c(surface, MeanNDSSI))

DeltaNDVIatMaxNDSSI

M <- DeltaNDVIatMaxNDSSI %>%
  ggplot( aes(x=MeanNDVI)) +
    geom_histogram( color="#e9ecef", alpha=0.6)

M

   
ggplot(DeltaNDVIatMaxNDSSI, aes(x = MeanNDVI)) + 
  geom_dotplot(binwidth = 0.1, dotsize = 0.25) + 
  scale_y_continuous(NULL, breaks = NULL)  + 
  xlim(0,1) +
  labs(x = "NDVI") +
  ggtitle("NDVI at month of maximum mean NDSSI")

And what if NDVI peak was one month sooner:

#extract NDVI value for each delta at one month earlier than max NDSSI value

DeltaNDVIatEarlyNDSSI <- DeltaMaxMin %>%
  select(Delta,MaxMeanNDSSImonth)

DeltaNDVIatEarlyNDSSI <- DeltaNDVIatEarlyNDSSI %>% 
  mutate(EarlyNDSSI = if_else(MaxMeanNDSSImonth == 1, 
                              12,
                              MaxMeanNDSSImonth-1)
         )

DeltaNDVIatEarlyNDSSI <- DeltaNDVIatEarlyNDSSI %>%
  select (-c(MaxMeanNDSSImonth))

DeltaMeansToJoin <- DeltaMeans %>%
  filter(surface == 'Land')

DeltaNDVIatEarlyNDSSI <- left_join(DeltaNDVIatEarlyNDSSI, DeltaMeansToJoin, 
                                 by = c('Delta', 'EarlyNDSSI' ='month'))
 
DeltaNDVIatEarlyNDSSI <- DeltaNDVIatEarlyNDSSI %>%
  select (-c(surface, MeanNDSSI))

DeltaNDVIatEarlyNDSSI

n <- DeltaNDVIatEarlyNDSSI %>%
  ggplot( aes(x=MeanNDVI)) +
    geom_histogram( color="#e9ecef", alpha=0.6)

n


ggplot(DeltaNDVIatEarlyNDSSI, aes(x = MeanNDVI)) + 
  geom_dotplot(binwidth = 0.1, dotsize = 0.25) + 
  scale_y_continuous(NULL, breaks = NULL)  + 
  xlim(0,1) +
  labs(x = "NDVI") +
  ggtitle("NDVI at 1 month earlier than maximum mean NDSSI")

How about the NDVI for months with the lowest NDSSI

#extract NDVI value for each delta a the month of min NDSSI value

DeltaNDVIatMinNDSSI <- DeltaMaxMin %>%
  select(Delta,MinMeanNDSSImonth)

DeltaMeansToJoin <- DeltaMeans %>%
  filter(surface == 'Land')

DeltaNDVIatMinNDSSI <- left_join(DeltaNDVIatMinNDSSI, DeltaMeansToJoin, 
                                 by = c('Delta', 'MinMeanNDSSImonth' ='month'))
 
DeltaNDVIatMinNDSSI <- DeltaNDVIatMinNDSSI %>%
  select (-c(surface, MeanNDSSI))

DeltaNDVIatMinNDSSI

p <- DeltaNDVIatMinNDSSI %>%
  ggplot( aes(x=MeanNDVI)) +
    geom_histogram( color="#e9ecef", alpha=0.6)

p


ggplot(DeltaNDVIatMinNDSSI, aes(x = MeanNDVI)) + 
  geom_dotplot(binwidth = 0.1, dotsize = 0.25) + 
  scale_y_continuous(NULL, breaks = NULL)  + 
  xlim(0,1) +
  labs(x = "NDVI") +
  ggtitle("NDVI at month of min mean NDSSI")

Just to explore the data a bit, here are the phase shifts/offsets against other measured parameters for each delta. The range, max and mean of NDVI and NDSSI is calculated from the timeseries, so it is really the max, min, and range of the monthly means (i.e., the maximum of the means, the minimum of the means, and the range of the mean). Offset is measured in months.

slice1 <- ggplot(DeltaMaxMin, aes(y = rangeNDVI, x = MaxOffset)) + geom_point() 
# + scale_x_discrete(limits = c(1:6), breaks = c(1:6)) + expand_limits(x = c(1,6)) 

slice2 <- ggplot(DeltaMaxMin, aes(y = rangeNDVI, x = rangeNDSSI)) + geom_point() 
slice3 <- ggplot(DeltaMaxMin, aes(y = rangeNDSSI, x = MaxOffset)) + geom_point() 

slice4 <- ggplot(DeltaMaxMin, aes(y = MaxMeanNDVI, x = rangeNDVI)) + geom_point() 
slice5 <- ggplot(DeltaMaxMin, aes(y = MaxMeanNDVI, x = rangeNDSSI)) + geom_point() 
slice6 <- ggplot(DeltaMaxMin, aes(y = MaxMeanNDVI, x = MaxOffset)) + geom_point() 

slice7 <- ggplot(DeltaMaxMin, aes(y = MaxMeanNDSSI, x = MaxMeanNDVI)) + geom_point() 
slice8 <- ggplot(DeltaMaxMin, aes(y = MaxMeanNDSSI, x = rangeNDSSI)) + geom_point() 
slice9 <- ggplot(DeltaMaxMin, aes(y = MaxMeanNDSSI, x = MaxOffset)) + geom_point() 
slice10 <- ggplot(DeltaMaxMin, aes(y = MaxMeanNDSSI, x = rangeNDVI)) + geom_point() 

grid.arrange(slice1, slice2, slice3, slice4, slice5, slice6, ncol=3)

grid.arrange(slice7, slice8, slice9, slice10, ncol=2)


#remove those panels
rm(slice1, slice2, slice3, slice4, slice5, slice6,slice7, slice8, slice9, slice10)

Join Latitude and longitude data

DeltaDatawLocations <- left_join(DeltaMaxMin, DeltaLocations, by = c("Delta" = "Deltas"))

DeltaDatawLocations <- DeltaDatawLocations %>%
  mutate(Absolute_Latitude= abs(Lat))

plot params vs lat

loc1 <- ggplot(DeltaDatawLocations, aes(y = Absolute_Latitude, x = MaxOffset)) + geom_point() 
loc2 <- ggplot(DeltaDatawLocations, aes(y = Absolute_Latitude, x = rangeNDSSI)) + geom_point() 
loc3 <- ggplot(DeltaDatawLocations, aes(y = Absolute_Latitude, x = rangeNDVI)) + geom_point() 
loc4 <- ggplot(DeltaDatawLocations, aes(y = Absolute_Latitude, x = MaxMeanNDVI)) + geom_point() 
loc5 <- ggplot(DeltaDatawLocations, aes(y = Absolute_Latitude, x = MaxMeanNDSSI)) + geom_point() 

grid.arrange(loc1, loc2, loc3, loc4, loc5, ncol=2)


loc1

#ggsave("loc1.pdf", width = 4, height = 4)
loc2

#ggsave("loc2.pdf", width = 4, height = 4)
loc3

#ggsave("loc3.pdf", width = 4, height = 4)

#remove those panels
rm(loc1, loc2, loc3, loc4, loc5)
#find the linear model 
DeltaOffset_lm <- lm( Absolute_Latitude ~ MaxOffset, data = DeltaDatawLocations) 

summary(DeltaOffset_lm)

Call:
lm(formula = Absolute_Latitude ~ MaxOffset, data = DeltaDatawLocations)

Residuals:
    Min      1Q  Median      3Q     Max 
-21.779  -8.983   0.775   8.377  20.772 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)   
(Intercept)   10.820      4.629   2.337  0.02653 * 
MaxOffset      4.079      1.157   3.525  0.00143 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 11.06 on 29 degrees of freedom
Multiple R-squared:  0.2999,    Adjusted R-squared:  0.2758 
F-statistic: 12.43 on 1 and 29 DF,  p-value: 0.001428
ggplot(DeltaDatawLocations, aes(x = Absolute_Latitude, y = MaxOffset)) + 
  geom_point() +
  geom_smooth(mapping = aes(x = Absolute_Latitude, y = MaxOffset, ), method=lm ) 

NA
NA
NA

Now for some maps of the data maps:

world <- ggplot() +
  borders("world", colour = "gray85", fill = "gray80") +
  theme_map() 

DeltaOffsetMap <- world +
  geom_point(aes(x = Lon, y = Lat, color = MaxOffset),
             data = DeltaDatawLocations, 
             size = 5) + scale_color_gradient( high = "red", low  = "yellow") +
  ggtitle("Offset Between NDVI peak on Land and NDSSI peak in water")

DeltaOffsetMap

#ggsave("DeltaOffsetMap.pdf", width = 6, height = 4)
world <- ggplot() +
  borders("world", colour = "gray85", fill = "gray80") +
  theme_map() 

DeltaNDVIrangeMap <- world +
  geom_point(aes(x = Lon, y = Lat, color = rangeNDVI),
             data = DeltaDatawLocations,
             size = 5) + scale_color_gradient( high = "red", low  = "yellow") + 
  ggtitle("NDVI range")


DeltaNDSSIrangeMap  <- world +
  geom_point(aes(x = Lon, y = Lat, color = rangeNDSSI),
             data = DeltaDatawLocations, 
             size = 5) + scale_color_gradient( high = "red", low = "yellow") + ggtitle("NDSSI range") 

DeltaNDVIrangeMap 

#ggsave("DeltaNDVIrangeMap.pdf", width = 6, height = 4)
DeltaNDSSIrangeMap

#ggsave("DeltaNDSSIrangeMap.pdf", width = 6, height = 4)

Let’s look at some of the timeseries To quantify the water, we use NDSSI. to quantify land, we use NDVI.

First here is the function to make the plots:

DeltaPlotter <- function(DeltaName) {
  #Counts each month
  numVeg <- DeltasCleaner %>%
    select(Delta, surface, month, ndvi) %>%
    filter(Delta == DeltaName & surface == "Land" & !is.na(ndvi)) %>%
    group_by(month) %>%
    summarize(n = n())
  
  numSed <- DeltasCleaner %>%
    select(Delta, surface, month, ndssi) %>%
    filter(Delta == DeltaName &
             surface == "Water" & !is.na(ndssi)) %>%
    group_by(month) %>%
    summarize(n = n())
  
  #Highlight the Maximum and Minimum Month for each delta, NDVI and NDSSI
  
  #LAND
  Veg <-
    ggplot(data = filter(DeltasCleaner, Delta == DeltaName &
                           surface == "Land")) +
    geom_boxplot(aes(x = month, y = ndvi, group = month)) +
    scale_x_discrete(limits = c(1:12), breaks = c(1:12)) +
    expand_limits(x = c(1, 12)) +
    ggtitle(DeltaName) +
    #geom_text(data = numVeg, aes(y = 1.05, x = month, label = n)) +
    geom_boxplot(
      data = filter(
        DeltasCleaner,
        Delta == DeltaName &
          surface == "Land" & month == DeltaMaxMin$MaxMeanNDVImonth[DeltaMaxMin$Delta == DeltaName] 
      ),
      aes(x = month, y = ndvi, group = month),
      fill = "green"
    ) +
    geom_boxplot(
      data = filter(
        DeltasCleaner,
        Delta == DeltaName & 
          surface == "Land" & month == DeltaMaxMin$MinMeanNDVImonth[DeltaMaxMin$Delta ==DeltaName]
      ),
      aes(x = month, y = ndvi, group = month),
      fill = "blue"
    )
  
  
  Sed <-
    ggplot(data = filter(DeltasCleaner, Delta == DeltaName &
                           surface == "Water")) +
    geom_boxplot(aes(x = month, y = ndssi, group = month)) +
    scale_x_discrete(limits = c(1:12), breaks = c(1:12)) +
    expand_limits(x = c(1, 12)) +
    #geom_text(data = numSed, aes(y = 1.05, x = month, label = n)) +
    geom_boxplot(
      data = filter(
        DeltasCleaner,
        Delta == DeltaName &
          surface == "Water" & month == DeltaMaxMin$MaxMeanNDSSImonth[DeltaMaxMin$Delta == DeltaName]
      ),
      aes(x = month, y = ndssi, group = month),
      fill = "green"
    ) +
    geom_boxplot(
      data = filter(
        DeltasCleaner,
        Delta == DeltaName &
          surface == "Water" & month == DeltaMaxMin$MinMeanNDSSImonth[DeltaMaxMin$Delta == DeltaName]
      ),
      aes(x = month, y = ndssi, group = month),
      fill = "blue"
    )
  
  return(grid.arrange(Veg, Sed, nrow = 2))
}

Here is are some examples:

DeltaPlotter("Parana")

DeltaPlotter("Magdalena")

DeltaPlotter("Ebro")

DeltaPlotter("Nile")

DeltaPlotter("Senegal")

DeltaPlotter("Orinoco")

DeltaPlotter("Godavari")

DeltaPlotter("Krishna")

Now for some work with GRDC discharge data:

#import the data (monthly means for 21 stations)
DeltasGRDC  <- read_csv("../data/GRDCstations.csv")
Parsed with column specification:
cols(
  Deltas = col_character(),
  GRDC_Station = col_double(),
  Time_Series_Length = col_character(),
  January = col_double(),
  February = col_double(),
  March = col_double(),
  April = col_double(),
  May = col_double(),
  June = col_double(),
  July = col_double(),
  August = col_double(),
  September = col_double(),
  October = col_double(),
  November = col_double(),
  December = col_double()
)
#calculate the mean of the monthly means
DeltasGRDCwMean <- DeltasGRDC %>% 
    rowwise() %>% 
    mutate(MMD=mean(c(January,February,March,April,
                    May,June,July,August,
                    September,October,November,December)))

DeltasGRDCwMean <- DeltasGRDCwMean %>% 
  rowwise() %>% 
  mutate(Range_Discharge = max(c(January,February,March,April,
              May,June,July,August,
              September,October,November,December)) - 
           min(c(January,February,March,April,
              May,June,July,August,
              September,October,November,December))
           )
        

#join to location data:
DeltawLocGRDC <- left_join(DeltaDatawLocations, DeltasGRDCwMean, by = c("Delta" = "Deltas"))


#plot mean of monthly means against NDSSI
ggplot(DeltawLocGRDC, aes(y = Range_Discharge, x = rangeNDSSI)) + geom_point() + scale_y_continuous(trans='log10')


#ggsave("GRDCNDSSI.pdf", width = 6, height = 4)

ggplot(DeltawLocGRDC, aes(y = Range_Discharge, x = MaxMeanNDSSI)) + geom_point() + scale_y_continuous(trans='log10')

#rename the months by numbers and tidy the GRDC data
DeltasDischarge <- DeltasGRDC %>%
  rename(Delta = Deltas,"1" = January, "2"= February, "3"= March, "4"= April,
         "5"=May, "6"=June, "7"=July, "8"= August, "9" = September, "10"=October, 
         "11"=November, "12"=December) %>%
  select(Delta, "1" , "2" , "3", "4","5", "6", "7", "8", "9", "10", "11", "12") %>%
  pivot_longer(-Delta, names_to = "month", values_to = "discharge")

#find max GRDC month for each delta
DeltaMaxDischarge <- 
  DeltasDischarge %>% 
  group_by(Delta) %>% 
  slice(which.max(discharge)) %>% 
  rename(MaxDischargeMonth = month, MaxDischarge = discharge) 


DeltaMaxDischarge$MaxDischargeMonth = as.numeric(DeltaMaxDischarge$MaxDischargeMonth)


#join with other delta data
DeltaMaxMinDischarge <- left_join(DeltaMaxMin, DeltaMaxDischarge, by = 'Delta')

#calculate offset
DeltaMaxMinDischarge <- DeltaMaxMinDischarge %>%
  mutate( DissOff = if_else(abs(MaxDischargeMonth - MaxMeanNDSSImonth) > 6,
                            12 - abs(MaxDischargeMonth - MaxMeanNDSSImonth),
                            abs(MaxDischargeMonth - MaxMeanNDSSImonth))
          )

#Compare offset with NDSSI (deltamaxmin$maxmeanNDSSImonth)
ggplot(DeltaMaxMinDischarge, aes(y = Delta, x = DissOff)) + geom_point() + 
  scale_x_discrete(limits = c(1:6), breaks = c(1:6)) +
  expand_limits(x = c(0,6))  + 
  ggtitle("DisOffset")

Look at GRDC data by latitude:

#join lat data
DeltaDatawLocations <- left_join(DeltaMaxMinDischarge, DeltaLocations, by = c("Delta" = "Deltas"))

DeltaDatawLocations <- DeltaDatawLocations %>%
  mutate(Absolute_Latitude= abs(Lat))

# plot offset on graph by lat
ggplot(DeltaDatawLocations, aes(x = Absolute_Latitude, y = DissOff)) + geom_point() +
  scale_color_gradient(low = "yellow", high = "red", na.value = NA) 

  #+ geom_smooth(mapping = aes(x = Absolute_Latitude, y = DissOff, ), method=lm ) 

#ggsave("DisOffset.pdf", width = 6, height = 4)

And now on a map:

#plot offset on map
DeltaDisOffsetMap <- world +
  geom_point(aes(x = Lon, y = Lat, color = DissOff),
             data = DeltaDatawLocations, 
             size = 5) + scale_color_gradient( high = "red", low  = "yellow") +
  ggtitle("Offset Between GRDC discharge peak and NDSSI peak in water")

DeltaDisOffsetMap

#ggsave("DeltaOffsetMap.pdf", width = 6, height = 4)
LS0tCnRpdGxlOiAiRGVsdGFzIE5vdGVib29rIgpvdXRwdXQ6CiAgaHRtbF9kb2N1bWVudDoKICAgIGRmX3ByaW50OiBwYWdlZAogIGh0bWxfbm90ZWJvb2s6IGRlZmF1bHQKICBwZGZfZG9jdW1lbnQ6IGRlZmF1bHQKICB3b3JkX2RvY3VtZW50OiBkZWZhdWx0CmVkaXRvcl9vcHRpb25zOiAKICBjaHVua19vdXRwdXRfdHlwZTogaW5saW5lCi0tLQoKTGV0J3MgaW1wb3J0IGF0IGFsbCBvZiB0aGUgZGF0YSB0aGF0IFNpbW9uIHB1bGxlZDoKCmBgYHtyfQpsaWJyYXJ5KHRpZHl2ZXJzZSkKbGlicmFyeShncmlkRXh0cmEpCmxpYnJhcnkobWFwcykKbGlicmFyeShnZ3RoZW1lcykKCkRlbHRhc0NsZWFuIDwtIHJlYWRfY3N2KCIuLi9kYXRhL291dC9kZWx0YXNfY2xlYW5fdjIuY3N2IikgCkRlbHRhTG9jYXRpb25zIDwtIHJlYWRfY3N2KCIuLi9kYXRhL0RlbHRhTG9jYXRpb25zLmNzdiIpCmBgYAoKQXMgYSByZW1pbmRlciwgZm9yIGVhY2ggb2YgdGhlIDQ3IGRlbHRhcyB0aGVyZSBhcmUgbWVhc3VyZW1lbnRzIG9mIExhbmQgJiBXYXRlciBhcmVhcyBhdCBVcHN0cmVhbSwgRG93bnN0cmVhbSBhbmQgJ01pZGRsZScgbG9jYXRpb25zIG9uIHRoZSBkZWx0YS4gV2UgZmlyc3QgbHVtcCBhbGwgdGhlIG9ic2VydmF0aW9ucyB0b2dldGhlciwgYW5kIGxvb2sgdG8gc2VlIHdoaWNoIERlbHRhcyBoYXZlIG1hbnkgb2JzZXJ2YXRpb25zOgoKYGBge3J9CiNjb3VudHMgcGVyIGRlbHRhCkRlbHRhQ291bnRzIDwtIGNvdW50KERlbHRhc0NsZWFuLCBEZWx0YSkKRGVsdGFDb3VudHMKYGBgCgoKTm93LCBieSBlYWNoIG1vbnRoLi4gd2hlcmUgdGhlIGNvbG9yYmFyIHJlcHJlc2VudHMgdGhlIG51bWJlciBvZiBvYnNlcnZhdGlvbnMgKG4pIGZvciBlYWNoIG1vbnRoIGZvciBhIGdpdmVuIGRlbHRhOgpgYGB7cn0KRGVsdGFPYnNQZXJNb250aCA8LSBjb3VudChEZWx0YXNDbGVhbiwgRGVsdGEsIG1vbnRoKQoKZ2dwbG90KERlbHRhT2JzUGVyTW9udGgsIGFlcyh5ID0gRGVsdGEsIHggPSBtb250aCwgZmlsbD1uKSkgKyBnZW9tX3RpbGUoKSArIAogIHNjYWxlX3hfZGlzY3JldGUobGltaXRzID0gYygxOjEyKSwgYnJlYWtzID0gYygxOjEyKSkgKwogIGV4cGFuZF9saW1pdHMoeCA9IGMoMSwxMikpICsgCiAgc2NhbGVfZmlsbF9ncmFkaWVudCggdHJhbnMgPSAnbG9nJyApCgpgYGAKCgpJbiB0aGUgYWJvdmUgaGVhdCBtYXAsIGRhcmsgY29sb3JzIChhbmQgbm8gY29sb3IpIHJlcHJlc2VudCBkYXRhIHBhdWNpdHkgKGFuZCBkYXRhIGdhcHMpLiBEZWx0YXMgd2l0aCBsaWdodCBjb2xvcnMgKGUuZy4sIHRoZSBQYXJhbmEsIE5pbGUsIEVicm8sIENvbG9yYWRvLCBCcmFobWFuaSkgaGF2ZSBsb3RzIG9mIGRhdGEsIHNwcmVhZCBvdXQgdGhyb3VnaCB0aGUgbW9udGhzIG9mIHRoZSB5ZWFyLgoKCkknbGwgcmVtb3ZlL3N1YnNldCB0aGUgZGVsdGFzIHdpdGggc3BhcnNlIGNvdmVyYWdlIChzcGVjaWZpY2FsbHksIG1vbnRocyB3aXRoIG5vIGNvdmVyYWdlKS4uLi4gCmBgYHtyfQoKIyBuZWVkIDEwIGRhdGEgcG9pbnRzIHBlciBtb250aCBmb3IgTkRTU0kgYW5kIE5EVkkKRW5vdWdoT2JzUGVyTW9udGggPC0gRGVsdGFzQ2xlYW4gJT4lIHVuZ3JvdXAoKSAlPiUKICBjb3VudChEZWx0YSwgbW9udGgsIHN1cmZhY2UpICU+JSAKICBncm91cF9ieShzdXJmYWNlKSAlPiUKICBmaWx0ZXIoIG4gPj0gNSkKCiNmaW5kIGRlbHRhcyBtaXNzaW5nIGEgZ2l2ZW4gbW9udGggb2Ygb2JzZXJ2YXRpb25zCkRlbHRhTW9udGhDb3VudHMgPC0gRW5vdWdoT2JzUGVyTW9udGggJT4lCiAgdW5ncm91cCgpICU+JQogIGNvdW50KERlbHRhKQoKIyBuZWVkIDEyIG1vbnRocyBvZiB3YXRlciBhbmQgbGFuZCBvYnMsIHNvIDI0IG1vIHRvdGFsCkVub3VnaE1vbnRocyA8LSBEZWx0YU1vbnRoQ291bnRzICU+JQogZmlsdGVyKCBuID09IDI0KQoKQ29tcGxldGVPYnNEZWx0YXMgPC0gcHVsbChFbm91Z2hNb250aHMsIERlbHRhKQoKI3JlbW92ZSB0aGVtCkRlbHRhc0NsZWFuZXIgPC0gRGVsdGFzQ2xlYW4gJT4lCiAgZmlsdGVyKERlbHRhICVpbiUgQ29tcGxldGVPYnNEZWx0YXMpCgojYWRkIHRoZSByZWFsIGRhdGVzIGluIG1vbnRoIGRhdGUgZm9ybWF0CkRlbHRhc0NsZWFuZXIkZGF0ZSA8LSBhcy5EYXRlKHBhc3RlKERlbHRhc0NsZWFuZXIkeWVhciwgRGVsdGFzQ2xlYW5lciRtb250aCwgIjAxIiwgc2VwPSItIiksICIlWS0lbS0lZCIpCgojcmVtb3ZlIGludGVybWVkaWF0ZSBkYXRhCnJtKENvbXBsZXRlT2JzRGVsdGFzLCBFbm91Z2hNb250aHMsIERlbHRhTW9udGhDb3VudHMpCgojRW5vdWdoTW9udGhzCmBgYAoKYW5kIGV4dHJhY3Qgc29tZSBtZXRyaWNzOyBzcGVjaWZpY2FsbHkgSSB3aWxsIG1ha2UgYSB0aW1lc2VyaWVzIG9mIE5EVkkgYW5kIE5EU1NJIGZvciBlYWNoIGRlbHRhIHVzaW5nIHRoZSBtZWFuIHZhbHVlIGZvciBlYWNoIG1vbnRoLgoKYGBge3IgIGluY2x1ZGUgPSBUUlVFfQojdGFrZSB0aGUgbWVhbiBORFZJIGFuZCBORFNTSSBmb3IgZWFjaCBtb250aCwgZm9yIGVhY2ggZGVsdGEKRGVsdGFNZWFucyA8LSBEZWx0YXNDbGVhbmVyICU+JQogIGdyb3VwX2J5KERlbHRhLCBtb250aCwgc3VyZmFjZSkgJT4lCiAgc3VtbWFyaXplKE1lYW5ORFZJID0gbWVhbihuZHZpLCBuYS5ybSA9IFRSVUUpLCBNZWFuTkRTU0kgPSBtZWFuKG5kc3NpLCBuYS5ybSA9IFRSVUUpKQoKI21ha2UgYSA5IGNvbHVtbiBkYXRhIGZyYW1lIHdpdGg6CiNkZWx0YSwgCiNtYXggYW5kIG1pbiBORFZJIG1vbnRoLCAKI05EU1NJIG1heCBhbmQgbWluIG1vbnRoLCAKI21heCBhbmQgbWluIHZhbHVlcyBmb3IgYm90aCBORFZJIGFuZCBORFNTSQoKIyMjIyMKCkRlbHRhTWF4TkRWSSA8LSAKICBEZWx0YU1lYW5zICU+JSAKICBmaWx0ZXIoc3VyZmFjZSA9PSAnTGFuZCcpICAlPiUgCiAgc2VsZWN0ICgtYyhNZWFuTkRTU0kpKSAlPiUgCiAgZ3JvdXBfYnkoRGVsdGEpICU+JSAKICBzbGljZSh3aGljaC5tYXgoTWVhbk5EVkkpKSAlPiUgCiAgcmVuYW1lKE1heE1lYW5ORFZJbW9udGggPSBtb250aCwgTWF4TWVhbk5EVkkgPSBNZWFuTkRWSSkKCkRlbHRhTWF4TkRTU0kgPC0gCiAgRGVsdGFNZWFucyAlPiUgCiAgZmlsdGVyKHN1cmZhY2UgPT0gJ1dhdGVyJykgICU+JSAKICBzZWxlY3QgKC1jKE1lYW5ORFZJKSkgJT4lIAogIGdyb3VwX2J5KERlbHRhKSAlPiUgCiAgc2xpY2Uod2hpY2gubWF4KE1lYW5ORFNTSSkpICU+JSAKICByZW5hbWUoTWF4TWVhbk5EU1NJbW9udGggPSBtb250aCwgTWF4TWVhbk5EU1NJID0gTWVhbk5EU1NJKQoKRGVsdGFNaW5ORFZJIDwtIAogIERlbHRhTWVhbnMgJT4lIAogIGZpbHRlcihzdXJmYWNlID09ICdMYW5kJykgICU+JSAKICBzZWxlY3QgKC1jKE1lYW5ORFNTSSkpICU+JSAKICBncm91cF9ieShEZWx0YSkgJT4lIAogIHNsaWNlKHdoaWNoLm1pbihNZWFuTkRWSSkpICU+JSAKICByZW5hbWUoTWluTWVhbk5EVkltb250aCA9IG1vbnRoLCBNaW5NZWFuTkRWSSA9IE1lYW5ORFZJKQoKRGVsdGFNaW5ORFNTSSA8LSAKICBEZWx0YU1lYW5zICU+JSAKICBmaWx0ZXIoc3VyZmFjZSA9PSAnV2F0ZXInKSAgJT4lIAogIHNlbGVjdCAoLWMoTWVhbk5EVkkpKSAlPiUgCiAgZ3JvdXBfYnkoRGVsdGEpICU+JSAKICBzbGljZSh3aGljaC5taW4oTWVhbk5EU1NJKSkgJT4lIAogIHJlbmFtZShNaW5NZWFuTkRTU0ltb250aCA9IG1vbnRoLCBNaW5NZWFuTkRTU0kgPSBNZWFuTkRTU0kpCgoKI2pvaW4gaW50byAxIGRhdGFmcmFtZQpEZWx0YU1heE1pbiA8LSBsZWZ0X2pvaW4oRGVsdGFNYXhORFZJLCBEZWx0YU1heE5EU1NJLCBieSA9ICdEZWx0YScpICU+JSAKICBsZWZ0X2pvaW4oLixEZWx0YU1pbk5EVkksIGJ5ID0gJ0RlbHRhJykgJT4lIAogIGxlZnRfam9pbiguLERlbHRhTWluTkRTU0ksIGJ5ID0gJ0RlbHRhJykKCiNyZW1vdmUgaW50ZXJtZWRpYXRlIGRhdGEKcm0oRGVsdGFNYXhORFZJLCBEZWx0YU1heE5EU1NJLCBEZWx0YU1pbk5EU1NJLERlbHRhTWluTkRWSSkKCmBgYAoKCkFuZCBub3cgd2UgY2FuIGxvb2sgYXQgdGhlIHBoYXNlIHNoaWZ0cyBiZXR3ZWVuIHRoZXNlIHR3byB0aW1lIHNlcmllcyAodGhlIHRpbWVzZXJpZXMgb2YgbWVhbiBORFZJIGFuZCBtZWFuIE5EU1NJKS4gSGVyZSBhcmUgdGhlIHBoYXNlIHNoaWZ0cyAoaW4gbW9udGgpIGZvciBlYWNoIGRlbHRhOgpgYGB7ciB9CiNjb21wYXJlIG9mZnNldApEZWx0YU1heE1pbiA8LSBtdXRhdGUoRGVsdGFNYXhNaW4sIAogICAgICAgICAgICAgICAgICAgICAgTWluT2Zmc2V0ID0gaWZfZWxzZShhYnMoTWluTWVhbk5EVkltb250aCAtIE1pbk1lYW5ORFNTSW1vbnRoKSA+IDYsIAogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAxMiAtIGFicyhNaW5NZWFuTkRWSW1vbnRoIC0gTWluTWVhbk5EU1NJbW9udGgpLAogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICBhYnMoTWluTWVhbk5EVkltb250aCAtIE1pbk1lYW5ORFNTSW1vbnRoKSksCiAgICAgICAgICAgICAgICAgICAgICBNYXhPZmZzZXQgPSBpZl9lbHNlKGFicyhNYXhNZWFuTkRWSW1vbnRoIC0gTWF4TWVhbk5EU1NJbW9udGgpID4gNiwgCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgMTIgLSBhYnMoTWF4TWVhbk5EVkltb250aCAtIE1heE1lYW5ORFNTSW1vbnRoKSwKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgYWJzKE1heE1lYW5ORFZJbW9udGggLSBNYXhNZWFuTkRTU0ltb250aCkpLAogICAgICAgICAgICAgICAgICAgICAgT2Zmc2V0RGlmZiA9IGFicyhNYXhPZmZzZXQgLSBNaW5PZmZzZXQpLAogICAgICAgICAgICAgICAgICAgICAgcmFuZ2VORFZJID0gKE1heE1lYW5ORFZJIC0gTWluTWVhbk5EVkkpLCAKICAgICAgICAgICAgICAgICAgICAgIHJhbmdlTkRTU0kgPSAoTWF4TWVhbk5EU1NJIC0gTWluTWVhbk5EU1NJKSwKICAgICAgICAgICAgICAgICAgICAgIGhhbGZQZXJpb2RORFZJID0gaWZfZWxzZShhYnMoTWF4TWVhbk5EVkltb250aCAtIE1pbk1lYW5ORFZJbW9udGgpID4gNiwgCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIDEyIC0gYWJzKE1heE1lYW5ORFZJbW9udGggLSBNaW5NZWFuTkRWSW1vbnRoKSwKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgYWJzKE1heE1lYW5ORFZJbW9udGggLSBNaW5NZWFuTkRWSW1vbnRoKSksCiAgICAgICAgICAgICAgICAgICAgICBoYWxmUGVyaW9kTkRTU0kgPSBpZl9lbHNlKGFicyhNYXhNZWFuTkRTU0ltb250aCAtIE1pbk1lYW5ORFNTSW1vbnRoKSA+IDYsIAogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAxMiAtIGFicyhNYXhNZWFuTkRTU0ltb250aCAtIE1pbk1lYW5ORFNTSW1vbnRoKSwKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgYWJzKE1heE1lYW5ORFNTSW1vbnRoIC0gTWluTWVhbk5EU1NJbW9udGgpKSwgKQoKIyBEZWx0YU1heE1pbiA8LSAKIyAgIERlbHRhTWF4TWluICAgJT4lCiMgICBzZWxlY3QgKGMoRGVsdGEsIE1pbk9mZnNldCwgTWF4T2Zmc2V0LCBPZmZzZXREaWZmLCByYW5nZU5EVkksIHJhbmdlTkRTU0ksTWF4TWVhbk5EU1NJLE1pbk1lYW5ORFNTSSxNYXhNZWFuTkRWSSxNaW5NZWFuTkRWSSkpIAoKRGVsdGFNYXhNaW4KCmdncGxvdChEZWx0YU1heE1pbiwgYWVzKHkgPSBEZWx0YSwgeCA9IE1heE9mZnNldCkpICsgZ2VvbV9wb2ludCgpICsgCiAgc2NhbGVfeF9kaXNjcmV0ZShsaW1pdHMgPSBjKDE6NiksIGJyZWFrcyA9IGMoMTo2KSkgKwogIGV4cGFuZF9saW1pdHMoeCA9IGMoMCw2KSkgICsgCiAgZ2d0aXRsZSgiTWF4T2Zmc2V0IikKCmdncGxvdChEZWx0YU1heE1pbiwgYWVzKHkgPSBEZWx0YSwgeCA9IE1pbk9mZnNldCkpICsgZ2VvbV9wb2ludCgpICsgCiAgc2NhbGVfeF9kaXNjcmV0ZShsaW1pdHMgPSBjKDE6NiksIGJyZWFrcyA9IGMoMTo2KSkgKwogIGV4cGFuZF9saW1pdHMoeCA9IGMoMCw2KSkgICsgCiAgZ2d0aXRsZSgiTWluT2Zmc2V0IikKCmdncGxvdChEZWx0YU1heE1pbiwgYWVzKHkgPSBEZWx0YSwgeCA9IE9mZnNldERpZmYpKSArIGdlb21fcG9pbnQoKSArIAogIHNjYWxlX3hfZGlzY3JldGUobGltaXRzID0gYygxOjYpLCBicmVha3MgPSBjKDE6NikpICsKICBleHBhbmRfbGltaXRzKHggPSBjKDAsNikpICArIAogIGdndGl0bGUoIk9mZnNldCBEaWZmZXJlbmNlIikKCmBgYAoKTm93IGxldCdzIGV4YW1pbmUgdGhlIGhpc3RvZ3JhbXMgb2YgYWxsIDMxIGRlbHRhcy4uLiBUaGUgbW9udGhzIHdpdGggdGhlIGdyZWF0ZXN0IG1lYW4gTkRWSSwgbW9udGhzIHdpdGggZ3JldGFlc3QgbWVhbiBORFNTSSwgdGhlIG1vbnRobHkgb2Zmc2V0LCBhbmQgdGhlIHNrZXcgb2YgdGhlIE5EU1NJIGFuZCBORFZJIHRpbWVzZXJpZXMuCgpgYGB7cn0KZ2dwbG90KERlbHRhTWF4TWluLCBhZXMoeCA9IE1heE1lYW5ORFZJbW9udGgpKSArIAogIGdlb21fZG90cGxvdChiaW53aWR0aCA9IDEsZG90c2l6ZSA9IDAuNSkgKyAKICBzY2FsZV95X2NvbnRpbnVvdXMoTlVMTCwgYnJlYWtzID0gTlVMTCkgKyAKICBzY2FsZV94X2Rpc2NyZXRlKGxpbWl0cyA9IGMoMToxMiksIGJyZWFrcyA9IGMoMToxMikpICsgCiAgbGFicyh4ID0gIk1vbnRoIikgKwogIGdndGl0bGUoIk1vbnRoIG9mIG1heGltdW0gbWVhbiBORFZJIikKCmdncGxvdChEZWx0YU1heE1pbiwgYWVzKHggPSBNYXhNZWFuTkRTU0ltb250aCkpICsgCiAgZ2VvbV9kb3RwbG90KGJpbndpZHRoID0gMSxkb3RzaXplID0gMC41KSArIAogIHNjYWxlX3lfY29udGludW91cyhOVUxMLCBicmVha3MgPSBOVUxMKSArIAogIHNjYWxlX3hfZGlzY3JldGUobGltaXRzID0gYygxOjEyKSwgYnJlYWtzID0gYygxOjEyKSkgKyAKICBsYWJzKHggPSAiTW9udGgiKSArCiAgZ2d0aXRsZSgiTW9udGggb2YgbWF4aW11bSBtZWFuIE5EU1NJIikKCmdncGxvdChEZWx0YU1heE1pbiwgYWVzKHggPSBNYXhPZmZzZXQpKSArIAogIGdlb21fZG90cGxvdChiaW53aWR0aCA9IDEsZG90c2l6ZSA9IDAuMjUpICsgCiAgc2NhbGVfeV9jb250aW51b3VzKE5VTEwsIGJyZWFrcyA9IE5VTEwpICsgCiAgc2NhbGVfeF9kaXNjcmV0ZShsaW1pdHMgPSBjKDA6NiksIGJyZWFrcyA9IGMoMDo2KSkgKyAKICBsYWJzKHggPSAiTW9udGhzIikgKwogIGdndGl0bGUoIk1vbnRocyBPZmZzZXQgYmV0d2VlbiBORFZJIGFuZCBORFNTSSIpCgpnZ3Bsb3QoRGVsdGFNYXhNaW4sIGFlcyh4ID0gaGFsZlBlcmlvZE5EVkkpKSArIAogIGdlb21fZG90cGxvdChiaW53aWR0aCA9IDEsZG90c2l6ZSA9IDAuMjUpICsgCiAgc2NhbGVfeV9jb250aW51b3VzKE5VTEwsIGJyZWFrcyA9IE5VTEwpICsgCiAgc2NhbGVfeF9kaXNjcmV0ZShsaW1pdHMgPSBjKDA6NiksIGJyZWFrcyA9IGMoMDo2KSkgKyAKICBsYWJzKHggPSAiTW9udGhzIikgKwogIGdndGl0bGUoImhhbGYgcGVyaW9kIE5EVkkgIikKCmdncGxvdChEZWx0YU1heE1pbiwgYWVzKHggPSBoYWxmUGVyaW9kTkRTU0kpKSArIAogIGdlb21fZG90cGxvdChiaW53aWR0aCA9IDEsZG90c2l6ZSA9IDAuMjUpICsgCiAgc2NhbGVfeV9jb250aW51b3VzKE5VTEwsIGJyZWFrcyA9IE5VTEwpICsgCiAgc2NhbGVfeF9kaXNjcmV0ZShsaW1pdHMgPSBjKDA6NiksIGJyZWFrcyA9IGMoMDo2KSkgKyAKICBsYWJzKHggPSAiTW9udGhzIikgKwogIGdndGl0bGUoImhhbGYgcGVyaW9kIE5EU1NJIikKCmBgYAoKT2ssIHNvIHRoZSBpZGVhIGlzIHRoYXQgcGVhayBORFNTSSBpcyBtb3JlIGVmZmVjdGl2ZSBpZiBpdCBvY2N1cnMgYXQgbW9kZXJhdGUgTkRWSSwgc28gbGV0J3MgbG9vayBhdCB0aGUgTkRWSSB2YWx1ZSBmb3IgdGhlIG1vbnRocyB3aXRoIHBlYWsgTkRTU0kuIApgYGB7cn0KI2V4dHJhY3QgTkRWSSB2YWx1ZSBmb3IgZWFjaCBkZWx0YSBhIHRoZSBtb250aCBvZiBtYXggTkRTU0kgdmFsdWUKCkRlbHRhTkRWSWF0TWF4TkRTU0kgPC0gRGVsdGFNYXhNaW4gJT4lCiAgc2VsZWN0KERlbHRhLE1heE1lYW5ORFNTSW1vbnRoKQoKRGVsdGFNZWFuc1RvSm9pbiA8LSBEZWx0YU1lYW5zICU+JQogIGZpbHRlcihzdXJmYWNlID09ICdMYW5kJykKCkRlbHRhTkRWSWF0TWF4TkRTU0kgPC0gbGVmdF9qb2luKERlbHRhTkRWSWF0TWF4TkRTU0ksIERlbHRhTWVhbnNUb0pvaW4sIAogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICBieSA9IGMoJ0RlbHRhJywgJ01heE1lYW5ORFNTSW1vbnRoJyA9J21vbnRoJykpCiAKRGVsdGFORFZJYXRNYXhORFNTSSA8LSBEZWx0YU5EVklhdE1heE5EU1NJICU+JQogIHNlbGVjdCAoLWMoc3VyZmFjZSwgTWVhbk5EU1NJKSkKCkRlbHRhTkRWSWF0TWF4TkRTU0kKCk0gPC0gRGVsdGFORFZJYXRNYXhORFNTSSAlPiUKICBnZ3Bsb3QoIGFlcyh4PU1lYW5ORFZJKSkgKwogICAgZ2VvbV9oaXN0b2dyYW0oIGNvbG9yPSIjZTllY2VmIiwgYWxwaGE9MC42KQoKTQogICAKZ2dwbG90KERlbHRhTkRWSWF0TWF4TkRTU0ksIGFlcyh4ID0gTWVhbk5EVkkpKSArIAogIGdlb21fZG90cGxvdChiaW53aWR0aCA9IDAuMSwgZG90c2l6ZSA9IDAuMjUpICsgCiAgc2NhbGVfeV9jb250aW51b3VzKE5VTEwsIGJyZWFrcyA9IE5VTEwpICArIAogIHhsaW0oMCwxKSArCiAgbGFicyh4ID0gIk5EVkkiKSArCiAgZ2d0aXRsZSgiTkRWSSBhdCBtb250aCBvZiBtYXhpbXVtIG1lYW4gTkRTU0kiKQpgYGAKCkFuZCB3aGF0IGlmIE5EVkkgcGVhayB3YXMgb25lIG1vbnRoIHNvb25lcjoKYGBge3J9CiNleHRyYWN0IE5EVkkgdmFsdWUgZm9yIGVhY2ggZGVsdGEgYXQgb25lIG1vbnRoIGVhcmxpZXIgdGhhbiBtYXggTkRTU0kgdmFsdWUKCkRlbHRhTkRWSWF0RWFybHlORFNTSSA8LSBEZWx0YU1heE1pbiAlPiUKICBzZWxlY3QoRGVsdGEsTWF4TWVhbk5EU1NJbW9udGgpCgpEZWx0YU5EVklhdEVhcmx5TkRTU0kgPC0gRGVsdGFORFZJYXRFYXJseU5EU1NJICU+JSAKICBtdXRhdGUoRWFybHlORFNTSSA9IGlmX2Vsc2UoTWF4TWVhbk5EU1NJbW9udGggPT0gMSwgCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIDEyLAogICAgICAgICAgICAgICAgICAgICAgICAgICAgICBNYXhNZWFuTkRTU0ltb250aC0xKQogICAgICAgICApCgpEZWx0YU5EVklhdEVhcmx5TkRTU0kgPC0gRGVsdGFORFZJYXRFYXJseU5EU1NJICU+JQogIHNlbGVjdCAoLWMoTWF4TWVhbk5EU1NJbW9udGgpKQoKRGVsdGFNZWFuc1RvSm9pbiA8LSBEZWx0YU1lYW5zICU+JQogIGZpbHRlcihzdXJmYWNlID09ICdMYW5kJykKCkRlbHRhTkRWSWF0RWFybHlORFNTSSA8LSBsZWZ0X2pvaW4oRGVsdGFORFZJYXRFYXJseU5EU1NJLCBEZWx0YU1lYW5zVG9Kb2luLCAKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgYnkgPSBjKCdEZWx0YScsICdFYXJseU5EU1NJJyA9J21vbnRoJykpCiAKRGVsdGFORFZJYXRFYXJseU5EU1NJIDwtIERlbHRhTkRWSWF0RWFybHlORFNTSSAlPiUKICBzZWxlY3QgKC1jKHN1cmZhY2UsIE1lYW5ORFNTSSkpCgpEZWx0YU5EVklhdEVhcmx5TkRTU0kKCm4gPC0gRGVsdGFORFZJYXRFYXJseU5EU1NJICU+JQogIGdncGxvdCggYWVzKHg9TWVhbk5EVkkpKSArCiAgICBnZW9tX2hpc3RvZ3JhbSggY29sb3I9IiNlOWVjZWYiLCBhbHBoYT0wLjYpCgpuCgpnZ3Bsb3QoRGVsdGFORFZJYXRFYXJseU5EU1NJLCBhZXMoeCA9IE1lYW5ORFZJKSkgKyAKICBnZW9tX2RvdHBsb3QoYmlud2lkdGggPSAwLjEsIGRvdHNpemUgPSAwLjI1KSArIAogIHNjYWxlX3lfY29udGludW91cyhOVUxMLCBicmVha3MgPSBOVUxMKSAgKyAKICB4bGltKDAsMSkgKwogIGxhYnMoeCA9ICJORFZJIikgKwogIGdndGl0bGUoIk5EVkkgYXQgMSBtb250aCBlYXJsaWVyIHRoYW4gbWF4aW11bSBtZWFuIE5EU1NJIikKYGBgCgoKSG93IGFib3V0IHRoZSBORFZJIGZvciBtb250aHMgd2l0aCB0aGUgbG93ZXN0IE5EU1NJCgpgYGB7cn0KI2V4dHJhY3QgTkRWSSB2YWx1ZSBmb3IgZWFjaCBkZWx0YSBhIHRoZSBtb250aCBvZiBtaW4gTkRTU0kgdmFsdWUKCkRlbHRhTkRWSWF0TWluTkRTU0kgPC0gRGVsdGFNYXhNaW4gJT4lCiAgc2VsZWN0KERlbHRhLE1pbk1lYW5ORFNTSW1vbnRoKQoKRGVsdGFNZWFuc1RvSm9pbiA8LSBEZWx0YU1lYW5zICU+JQogIGZpbHRlcihzdXJmYWNlID09ICdMYW5kJykKCkRlbHRhTkRWSWF0TWluTkRTU0kgPC0gbGVmdF9qb2luKERlbHRhTkRWSWF0TWluTkRTU0ksIERlbHRhTWVhbnNUb0pvaW4sIAogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICBieSA9IGMoJ0RlbHRhJywgJ01pbk1lYW5ORFNTSW1vbnRoJyA9J21vbnRoJykpCiAKRGVsdGFORFZJYXRNaW5ORFNTSSA8LSBEZWx0YU5EVklhdE1pbk5EU1NJICU+JQogIHNlbGVjdCAoLWMoc3VyZmFjZSwgTWVhbk5EU1NJKSkKCkRlbHRhTkRWSWF0TWluTkRTU0kKCnAgPC0gRGVsdGFORFZJYXRNaW5ORFNTSSAlPiUKICBnZ3Bsb3QoIGFlcyh4PU1lYW5ORFZJKSkgKwogICAgZ2VvbV9oaXN0b2dyYW0oIGNvbG9yPSIjZTllY2VmIiwgYWxwaGE9MC42KQoKcAoKZ2dwbG90KERlbHRhTkRWSWF0TWluTkRTU0ksIGFlcyh4ID0gTWVhbk5EVkkpKSArIAogIGdlb21fZG90cGxvdChiaW53aWR0aCA9IDAuMSwgZG90c2l6ZSA9IDAuMjUpICsgCiAgc2NhbGVfeV9jb250aW51b3VzKE5VTEwsIGJyZWFrcyA9IE5VTEwpICArIAogIHhsaW0oMCwxKSArCiAgbGFicyh4ID0gIk5EVkkiKSArCiAgZ2d0aXRsZSgiTkRWSSBhdCBtb250aCBvZiBtaW4gbWVhbiBORFNTSSIpCmBgYAoKCkp1c3QgdG8gZXhwbG9yZSB0aGUgZGF0YSBhIGJpdCwgaGVyZSBhcmUgdGhlIHBoYXNlIHNoaWZ0cy9vZmZzZXRzIGFnYWluc3Qgb3RoZXIgbWVhc3VyZWQgcGFyYW1ldGVycyBmb3IgZWFjaCBkZWx0YS4gVGhlIHJhbmdlLCBtYXggYW5kIG1lYW4gb2YgTkRWSSBhbmQgTkRTU0kgaXMgY2FsY3VsYXRlZCBmcm9tIHRoZSB0aW1lc2VyaWVzLCBzbyBpdCBpcyByZWFsbHkgdGhlIG1heCwgbWluLCBhbmQgcmFuZ2Ugb2YgdGhlIG1vbnRobHkgbWVhbnMgKGkuZS4sIHRoZSBtYXhpbXVtIG9mIHRoZSBtZWFucywgdGhlIG1pbmltdW0gb2YgdGhlIG1lYW5zLCBhbmQgdGhlIHJhbmdlIG9mIHRoZSBtZWFuKS4gT2Zmc2V0IGlzIG1lYXN1cmVkIGluIG1vbnRocy4KCmBgYHtyfQpzbGljZTEgPC0gZ2dwbG90KERlbHRhTWF4TWluLCBhZXMoeSA9IHJhbmdlTkRWSSwgeCA9IE1heE9mZnNldCkpICsgZ2VvbV9wb2ludCgpIAojICsgc2NhbGVfeF9kaXNjcmV0ZShsaW1pdHMgPSBjKDE6NiksIGJyZWFrcyA9IGMoMTo2KSkgKyBleHBhbmRfbGltaXRzKHggPSBjKDEsNikpIAoKc2xpY2UyIDwtIGdncGxvdChEZWx0YU1heE1pbiwgYWVzKHkgPSByYW5nZU5EVkksIHggPSByYW5nZU5EU1NJKSkgKyBnZW9tX3BvaW50KCkgCnNsaWNlMyA8LSBnZ3Bsb3QoRGVsdGFNYXhNaW4sIGFlcyh5ID0gcmFuZ2VORFNTSSwgeCA9IE1heE9mZnNldCkpICsgZ2VvbV9wb2ludCgpIAoKc2xpY2U0IDwtIGdncGxvdChEZWx0YU1heE1pbiwgYWVzKHkgPSBNYXhNZWFuTkRWSSwgeCA9IHJhbmdlTkRWSSkpICsgZ2VvbV9wb2ludCgpIApzbGljZTUgPC0gZ2dwbG90KERlbHRhTWF4TWluLCBhZXMoeSA9IE1heE1lYW5ORFZJLCB4ID0gcmFuZ2VORFNTSSkpICsgZ2VvbV9wb2ludCgpIApzbGljZTYgPC0gZ2dwbG90KERlbHRhTWF4TWluLCBhZXMoeSA9IE1heE1lYW5ORFZJLCB4ID0gTWF4T2Zmc2V0KSkgKyBnZW9tX3BvaW50KCkgCgpzbGljZTcgPC0gZ2dwbG90KERlbHRhTWF4TWluLCBhZXMoeSA9IE1heE1lYW5ORFNTSSwgeCA9IE1heE1lYW5ORFZJKSkgKyBnZW9tX3BvaW50KCkgCnNsaWNlOCA8LSBnZ3Bsb3QoRGVsdGFNYXhNaW4sIGFlcyh5ID0gTWF4TWVhbk5EU1NJLCB4ID0gcmFuZ2VORFNTSSkpICsgZ2VvbV9wb2ludCgpIApzbGljZTkgPC0gZ2dwbG90KERlbHRhTWF4TWluLCBhZXMoeSA9IE1heE1lYW5ORFNTSSwgeCA9IE1heE9mZnNldCkpICsgZ2VvbV9wb2ludCgpIApzbGljZTEwIDwtIGdncGxvdChEZWx0YU1heE1pbiwgYWVzKHkgPSBNYXhNZWFuTkRTU0ksIHggPSByYW5nZU5EVkkpKSArIGdlb21fcG9pbnQoKSAKCmdyaWQuYXJyYW5nZShzbGljZTEsIHNsaWNlMiwgc2xpY2UzLCBzbGljZTQsIHNsaWNlNSwgc2xpY2U2LCBuY29sPTMpCmdyaWQuYXJyYW5nZShzbGljZTcsIHNsaWNlOCwgc2xpY2U5LCBzbGljZTEwLCBuY29sPTIpCgojcmVtb3ZlIHRob3NlIHBhbmVscwpybShzbGljZTEsIHNsaWNlMiwgc2xpY2UzLCBzbGljZTQsIHNsaWNlNSwgc2xpY2U2LHNsaWNlNywgc2xpY2U4LCBzbGljZTksIHNsaWNlMTApCgpgYGAKCgpKb2luIExhdGl0dWRlIGFuZCBsb25naXR1ZGUgZGF0YQoKYGBge3J9CkRlbHRhRGF0YXdMb2NhdGlvbnMgPC0gbGVmdF9qb2luKERlbHRhTWF4TWluLCBEZWx0YUxvY2F0aW9ucywgYnkgPSBjKCJEZWx0YSIgPSAiRGVsdGFzIikpCgpEZWx0YURhdGF3TG9jYXRpb25zIDwtIERlbHRhRGF0YXdMb2NhdGlvbnMgJT4lCiAgbXV0YXRlKEFic29sdXRlX0xhdGl0dWRlPSBhYnMoTGF0KSkKYGBgCgpwbG90IHBhcmFtcyB2cyBsYXQKYGBge3J9CmxvYzEgPC0gZ2dwbG90KERlbHRhRGF0YXdMb2NhdGlvbnMsIGFlcyh5ID0gQWJzb2x1dGVfTGF0aXR1ZGUsIHggPSBNYXhPZmZzZXQpKSArIGdlb21fcG9pbnQoKSAKbG9jMiA8LSBnZ3Bsb3QoRGVsdGFEYXRhd0xvY2F0aW9ucywgYWVzKHkgPSBBYnNvbHV0ZV9MYXRpdHVkZSwgeCA9IHJhbmdlTkRTU0kpKSArIGdlb21fcG9pbnQoKSAKbG9jMyA8LSBnZ3Bsb3QoRGVsdGFEYXRhd0xvY2F0aW9ucywgYWVzKHkgPSBBYnNvbHV0ZV9MYXRpdHVkZSwgeCA9IHJhbmdlTkRWSSkpICsgZ2VvbV9wb2ludCgpIApsb2M0IDwtIGdncGxvdChEZWx0YURhdGF3TG9jYXRpb25zLCBhZXMoeSA9IEFic29sdXRlX0xhdGl0dWRlLCB4ID0gTWF4TWVhbk5EVkkpKSArIGdlb21fcG9pbnQoKSAKbG9jNSA8LSBnZ3Bsb3QoRGVsdGFEYXRhd0xvY2F0aW9ucywgYWVzKHkgPSBBYnNvbHV0ZV9MYXRpdHVkZSwgeCA9IE1heE1lYW5ORFNTSSkpICsgZ2VvbV9wb2ludCgpIAoKZ3JpZC5hcnJhbmdlKGxvYzEsIGxvYzIsIGxvYzMsIGxvYzQsIGxvYzUsIG5jb2w9MikKCmxvYzEKI2dnc2F2ZSgibG9jMS5wZGYiLCB3aWR0aCA9IDQsIGhlaWdodCA9IDQpCmxvYzIKI2dnc2F2ZSgibG9jMi5wZGYiLCB3aWR0aCA9IDQsIGhlaWdodCA9IDQpCmxvYzMKI2dnc2F2ZSgibG9jMy5wZGYiLCB3aWR0aCA9IDQsIGhlaWdodCA9IDQpCgojcmVtb3ZlIHRob3NlIHBhbmVscwpybShsb2MxLCBsb2MyLCBsb2MzLCBsb2M0LCBsb2M1KQpgYGAKCmBgYHtyfQojZmluZCB0aGUgbGluZWFyIG1vZGVsIApEZWx0YU9mZnNldF9sbSA8LSBsbSggQWJzb2x1dGVfTGF0aXR1ZGUgfiBNYXhPZmZzZXQsIGRhdGEgPSBEZWx0YURhdGF3TG9jYXRpb25zKSAKCnN1bW1hcnkoRGVsdGFPZmZzZXRfbG0pCgpnZ3Bsb3QoRGVsdGFEYXRhd0xvY2F0aW9ucywgYWVzKHggPSBBYnNvbHV0ZV9MYXRpdHVkZSwgeSA9IE1heE9mZnNldCkpICsgCiAgZ2VvbV9wb2ludCgpICsKICBnZW9tX3Ntb290aChtYXBwaW5nID0gYWVzKHggPSBBYnNvbHV0ZV9MYXRpdHVkZSwgeSA9IE1heE9mZnNldCwgKSwgbWV0aG9kPWxtICkgCgoKCmBgYAoKCk5vdyBmb3Igc29tZSBtYXBzIG9mIHRoZSBkYXRhIG1hcHM6CgpgYGB7cn0Kd29ybGQgPC0gZ2dwbG90KCkgKwogIGJvcmRlcnMoIndvcmxkIiwgY29sb3VyID0gImdyYXk4NSIsIGZpbGwgPSAiZ3JheTgwIikgKwogIHRoZW1lX21hcCgpIAoKRGVsdGFPZmZzZXRNYXAgPC0gd29ybGQgKwogIGdlb21fcG9pbnQoYWVzKHggPSBMb24sIHkgPSBMYXQsIGNvbG9yID0gTWF4T2Zmc2V0KSwKICAgICAgICAgICAgIGRhdGEgPSBEZWx0YURhdGF3TG9jYXRpb25zLCAKICAgICAgICAgICAgIHNpemUgPSA1KSArIHNjYWxlX2NvbG9yX2dyYWRpZW50KCBoaWdoID0gInJlZCIsIGxvdyAgPSAieWVsbG93IikgKwogIGdndGl0bGUoIk9mZnNldCBCZXR3ZWVuIE5EVkkgcGVhayBvbiBMYW5kIGFuZCBORFNTSSBwZWFrIGluIHdhdGVyIikKCkRlbHRhT2Zmc2V0TWFwCiNnZ3NhdmUoIkRlbHRhT2Zmc2V0TWFwLnBkZiIsIHdpZHRoID0gNiwgaGVpZ2h0ID0gNCkKCmBgYAoKCmBgYHtyfQp3b3JsZCA8LSBnZ3Bsb3QoKSArCiAgYm9yZGVycygid29ybGQiLCBjb2xvdXIgPSAiZ3JheTg1IiwgZmlsbCA9ICJncmF5ODAiKSArCiAgdGhlbWVfbWFwKCkgCgpEZWx0YU5EVklyYW5nZU1hcCA8LSB3b3JsZCArCiAgZ2VvbV9wb2ludChhZXMoeCA9IExvbiwgeSA9IExhdCwgY29sb3IgPSByYW5nZU5EVkkpLAogICAgICAgICAgICAgZGF0YSA9IERlbHRhRGF0YXdMb2NhdGlvbnMsCiAgICAgICAgICAgICBzaXplID0gNSkgKyBzY2FsZV9jb2xvcl9ncmFkaWVudCggaGlnaCA9ICJyZWQiLCBsb3cgID0gInllbGxvdyIpICsgCiAgZ2d0aXRsZSgiTkRWSSByYW5nZSIpCgoKRGVsdGFORFNTSXJhbmdlTWFwICA8LSB3b3JsZCArCiAgZ2VvbV9wb2ludChhZXMoeCA9IExvbiwgeSA9IExhdCwgY29sb3IgPSByYW5nZU5EU1NJKSwKICAgICAgICAgICAgIGRhdGEgPSBEZWx0YURhdGF3TG9jYXRpb25zLCAKICAgICAgICAgICAgIHNpemUgPSA1KSArIHNjYWxlX2NvbG9yX2dyYWRpZW50KCBoaWdoID0gInJlZCIsIGxvdyA9ICJ5ZWxsb3ciKSArIGdndGl0bGUoIk5EU1NJIHJhbmdlIikgCgpEZWx0YU5EVklyYW5nZU1hcCAKI2dnc2F2ZSgiRGVsdGFORFZJcmFuZ2VNYXAucGRmIiwgd2lkdGggPSA2LCBoZWlnaHQgPSA0KQpEZWx0YU5EU1NJcmFuZ2VNYXAKI2dnc2F2ZSgiRGVsdGFORFNTSXJhbmdlTWFwLnBkZiIsIHdpZHRoID0gNiwgaGVpZ2h0ID0gNCkKCmBgYAoKTGV0J3MgbG9vayBhdCBzb21lIG9mIHRoZSB0aW1lc2VyaWVzIApUbyBxdWFudGlmeSB0aGUgd2F0ZXIsIHdlIHVzZSBORFNTSS4gdG8gcXVhbnRpZnkgbGFuZCwgd2UgdXNlIE5EVkkuCgpGaXJzdCBoZXJlIGlzIHRoZSBmdW5jdGlvbiB0byBtYWtlIHRoZSBwbG90czoKCmBgYHtyfQpEZWx0YVBsb3R0ZXIgPC0gZnVuY3Rpb24oRGVsdGFOYW1lKSB7CiAgI0NvdW50cyBlYWNoIG1vbnRoCiAgbnVtVmVnIDwtIERlbHRhc0NsZWFuZXIgJT4lCiAgICBzZWxlY3QoRGVsdGEsIHN1cmZhY2UsIG1vbnRoLCBuZHZpKSAlPiUKICAgIGZpbHRlcihEZWx0YSA9PSBEZWx0YU5hbWUgJiBzdXJmYWNlID09ICJMYW5kIiAmICFpcy5uYShuZHZpKSkgJT4lCiAgICBncm91cF9ieShtb250aCkgJT4lCiAgICBzdW1tYXJpemUobiA9IG4oKSkKICAKICBudW1TZWQgPC0gRGVsdGFzQ2xlYW5lciAlPiUKICAgIHNlbGVjdChEZWx0YSwgc3VyZmFjZSwgbW9udGgsIG5kc3NpKSAlPiUKICAgIGZpbHRlcihEZWx0YSA9PSBEZWx0YU5hbWUgJgogICAgICAgICAgICAgc3VyZmFjZSA9PSAiV2F0ZXIiICYgIWlzLm5hKG5kc3NpKSkgJT4lCiAgICBncm91cF9ieShtb250aCkgJT4lCiAgICBzdW1tYXJpemUobiA9IG4oKSkKICAKICAjSGlnaGxpZ2h0IHRoZSBNYXhpbXVtIGFuZCBNaW5pbXVtIE1vbnRoIGZvciBlYWNoIGRlbHRhLCBORFZJIGFuZCBORFNTSQogIAogICNMQU5ECiAgVmVnIDwtCiAgICBnZ3Bsb3QoZGF0YSA9IGZpbHRlcihEZWx0YXNDbGVhbmVyLCBEZWx0YSA9PSBEZWx0YU5hbWUgJgogICAgICAgICAgICAgICAgICAgICAgICAgICBzdXJmYWNlID09ICJMYW5kIikpICsKICAgIGdlb21fYm94cGxvdChhZXMoeCA9IG1vbnRoLCB5ID0gbmR2aSwgZ3JvdXAgPSBtb250aCkpICsKICAgIHNjYWxlX3hfZGlzY3JldGUobGltaXRzID0gYygxOjEyKSwgYnJlYWtzID0gYygxOjEyKSkgKwogICAgZXhwYW5kX2xpbWl0cyh4ID0gYygxLCAxMikpICsKICAgIGdndGl0bGUoRGVsdGFOYW1lKSArCiAgICAjZ2VvbV90ZXh0KGRhdGEgPSBudW1WZWcsIGFlcyh5ID0gMS4wNSwgeCA9IG1vbnRoLCBsYWJlbCA9IG4pKSArCiAgICBnZW9tX2JveHBsb3QoCiAgICAgIGRhdGEgPSBmaWx0ZXIoCiAgICAgICAgRGVsdGFzQ2xlYW5lciwKICAgICAgICBEZWx0YSA9PSBEZWx0YU5hbWUgJgogICAgICAgICAgc3VyZmFjZSA9PSAiTGFuZCIgJiBtb250aCA9PSBEZWx0YU1heE1pbiRNYXhNZWFuTkRWSW1vbnRoW0RlbHRhTWF4TWluJERlbHRhID09IERlbHRhTmFtZV0gCiAgICAgICksCiAgICAgIGFlcyh4ID0gbW9udGgsIHkgPSBuZHZpLCBncm91cCA9IG1vbnRoKSwKICAgICAgZmlsbCA9ICJncmVlbiIKICAgICkgKwogICAgZ2VvbV9ib3hwbG90KAogICAgICBkYXRhID0gZmlsdGVyKAogICAgICAgIERlbHRhc0NsZWFuZXIsCiAgICAgICAgRGVsdGEgPT0gRGVsdGFOYW1lICYgCiAgICAgICAgICBzdXJmYWNlID09ICJMYW5kIiAmIG1vbnRoID09IERlbHRhTWF4TWluJE1pbk1lYW5ORFZJbW9udGhbRGVsdGFNYXhNaW4kRGVsdGEgPT1EZWx0YU5hbWVdCiAgICAgICksCiAgICAgIGFlcyh4ID0gbW9udGgsIHkgPSBuZHZpLCBncm91cCA9IG1vbnRoKSwKICAgICAgZmlsbCA9ICJibHVlIgogICAgKQogIAogIAogIFNlZCA8LQogICAgZ2dwbG90KGRhdGEgPSBmaWx0ZXIoRGVsdGFzQ2xlYW5lciwgRGVsdGEgPT0gRGVsdGFOYW1lICYKICAgICAgICAgICAgICAgICAgICAgICAgICAgc3VyZmFjZSA9PSAiV2F0ZXIiKSkgKwogICAgZ2VvbV9ib3hwbG90KGFlcyh4ID0gbW9udGgsIHkgPSBuZHNzaSwgZ3JvdXAgPSBtb250aCkpICsKICAgIHNjYWxlX3hfZGlzY3JldGUobGltaXRzID0gYygxOjEyKSwgYnJlYWtzID0gYygxOjEyKSkgKwogICAgZXhwYW5kX2xpbWl0cyh4ID0gYygxLCAxMikpICsKICAgICNnZW9tX3RleHQoZGF0YSA9IG51bVNlZCwgYWVzKHkgPSAxLjA1LCB4ID0gbW9udGgsIGxhYmVsID0gbikpICsKICAgIGdlb21fYm94cGxvdCgKICAgICAgZGF0YSA9IGZpbHRlcigKICAgICAgICBEZWx0YXNDbGVhbmVyLAogICAgICAgIERlbHRhID09IERlbHRhTmFtZSAmCiAgICAgICAgICBzdXJmYWNlID09ICJXYXRlciIgJiBtb250aCA9PSBEZWx0YU1heE1pbiRNYXhNZWFuTkRTU0ltb250aFtEZWx0YU1heE1pbiREZWx0YSA9PSBEZWx0YU5hbWVdCiAgICAgICksCiAgICAgIGFlcyh4ID0gbW9udGgsIHkgPSBuZHNzaSwgZ3JvdXAgPSBtb250aCksCiAgICAgIGZpbGwgPSAiZ3JlZW4iCiAgICApICsKICAgIGdlb21fYm94cGxvdCgKICAgICAgZGF0YSA9IGZpbHRlcigKICAgICAgICBEZWx0YXNDbGVhbmVyLAogICAgICAgIERlbHRhID09IERlbHRhTmFtZSAmCiAgICAgICAgICBzdXJmYWNlID09ICJXYXRlciIgJiBtb250aCA9PSBEZWx0YU1heE1pbiRNaW5NZWFuTkRTU0ltb250aFtEZWx0YU1heE1pbiREZWx0YSA9PSBEZWx0YU5hbWVdCiAgICAgICksCiAgICAgIGFlcyh4ID0gbW9udGgsIHkgPSBuZHNzaSwgZ3JvdXAgPSBtb250aCksCiAgICAgIGZpbGwgPSAiYmx1ZSIKICAgICkKICAKICByZXR1cm4oZ3JpZC5hcnJhbmdlKFZlZywgU2VkLCBucm93ID0gMikpCn0KYGBgCgoKSGVyZSBpcyBhcmUgc29tZSBleGFtcGxlczoKCiogVGhlIHBlYWtzIGluIGJvdGggdGltZXNlcmllcyBzaGlmdCBhcm91bmQgZGVwZW5kaW5nIG9uIHRoZSBkZWx0YToKICAgKyBsb29rIGF0IHRoZSBjb3JyZWxhdGlvbiBpbiB0aGUgT3Jpbm9jbyBuYWQgU2VuZWdhbAogICArIFRoZSBhbnRpY29ycmVsYXRpb24gaW4gdGhlIFBhcmFuYSBhbmQgRWJybywKICAgKyBUaGUgc2xpZ2h0IHBoYXNlIHNoaWZ0IGluIHRoZSBNYWdkYWxlbmEuCgpgYGB7cn0KRGVsdGFQbG90dGVyKCJQYXJhbmEiKQpEZWx0YVBsb3R0ZXIoIk1hZ2RhbGVuYSIpCkRlbHRhUGxvdHRlcigiRWJybyIpCkRlbHRhUGxvdHRlcigiTmlsZSIpCkRlbHRhUGxvdHRlcigiU2VuZWdhbCIpCkRlbHRhUGxvdHRlcigiT3Jpbm9jbyIpCkRlbHRhUGxvdHRlcigiR29kYXZhcmkiKQpEZWx0YVBsb3R0ZXIoIktyaXNobmEiKQpgYGAKCgpOb3cgZm9yIHNvbWUgd29yayB3aXRoIEdSREMgZGlzY2hhcmdlIGRhdGE6CmBgYHtyfQojaW1wb3J0IHRoZSBkYXRhIChtb250aGx5IG1lYW5zIGZvciAyMSBzdGF0aW9ucykKRGVsdGFzR1JEQyAgPC0gcmVhZF9jc3YoIi4uL2RhdGEvR1JEQ3N0YXRpb25zLmNzdiIpCgojY2FsY3VsYXRlIHRoZSBtZWFuIG9mIHRoZSBtb250aGx5IG1lYW5zCkRlbHRhc0dSREN3TWVhbiA8LSBEZWx0YXNHUkRDICU+JSAKICAgIHJvd3dpc2UoKSAlPiUgCiAgICBtdXRhdGUoTU1EPW1lYW4oYyhKYW51YXJ5LEZlYnJ1YXJ5LE1hcmNoLEFwcmlsLAogICAgICAgICAgICAgICAgICAgIE1heSxKdW5lLEp1bHksQXVndXN0LAogICAgICAgICAgICAgICAgICAgIFNlcHRlbWJlcixPY3RvYmVyLE5vdmVtYmVyLERlY2VtYmVyKSkpCgpEZWx0YXNHUkRDd01lYW4gPC0gRGVsdGFzR1JEQ3dNZWFuICU+JSAKICByb3d3aXNlKCkgJT4lIAogIG11dGF0ZShSYW5nZV9EaXNjaGFyZ2UgPSBtYXgoYyhKYW51YXJ5LEZlYnJ1YXJ5LE1hcmNoLEFwcmlsLAogICAgICAgICAgICAgIE1heSxKdW5lLEp1bHksQXVndXN0LAogICAgICAgICAgICAgIFNlcHRlbWJlcixPY3RvYmVyLE5vdmVtYmVyLERlY2VtYmVyKSkgLSAKICAgICAgICAgICBtaW4oYyhKYW51YXJ5LEZlYnJ1YXJ5LE1hcmNoLEFwcmlsLAogICAgICAgICAgICAgIE1heSxKdW5lLEp1bHksQXVndXN0LAogICAgICAgICAgICAgIFNlcHRlbWJlcixPY3RvYmVyLE5vdmVtYmVyLERlY2VtYmVyKSkKICAgICAgICAgICApCiAgICAgICAgCgojam9pbiB0byBsb2NhdGlvbiBkYXRhOgpEZWx0YXdMb2NHUkRDIDwtIGxlZnRfam9pbihEZWx0YURhdGF3TG9jYXRpb25zLCBEZWx0YXNHUkRDd01lYW4sIGJ5ID0gYygiRGVsdGEiID0gIkRlbHRhcyIpKQoKCiNwbG90IG1lYW4gb2YgbW9udGhseSBtZWFucyBhZ2FpbnN0IE5EU1NJCmdncGxvdChEZWx0YXdMb2NHUkRDLCBhZXMoeSA9IFJhbmdlX0Rpc2NoYXJnZSwgeCA9IHJhbmdlTkRTU0kpKSArIGdlb21fcG9pbnQoKSArIHNjYWxlX3lfY29udGludW91cyh0cmFucz0nbG9nMTAnKQoKI2dnc2F2ZSgiR1JEQ05EU1NJLnBkZiIsIHdpZHRoID0gNiwgaGVpZ2h0ID0gNCkKCmdncGxvdChEZWx0YXdMb2NHUkRDLCBhZXMoeSA9IFJhbmdlX0Rpc2NoYXJnZSwgeCA9IE1heE1lYW5ORFNTSSkpICsgZ2VvbV9wb2ludCgpICsgc2NhbGVfeV9jb250aW51b3VzKHRyYW5zPSdsb2cxMCcpCmBgYAoKCmBgYHtyfQojcmVuYW1lIHRoZSBtb250aHMgYnkgbnVtYmVycyBhbmQgdGlkeSB0aGUgR1JEQyBkYXRhCkRlbHRhc0Rpc2NoYXJnZSA8LSBEZWx0YXNHUkRDICU+JQogIHJlbmFtZShEZWx0YSA9IERlbHRhcywiMSIgPSBKYW51YXJ5LCAiMiI9IEZlYnJ1YXJ5LCAiMyI9IE1hcmNoLCAiNCI9IEFwcmlsLAogICAgICAgICAiNSI9TWF5LCAiNiI9SnVuZSwgIjciPUp1bHksICI4Ij0gQXVndXN0LCAiOSIgPSBTZXB0ZW1iZXIsICIxMCI9T2N0b2JlciwgCiAgICAgICAgICIxMSI9Tm92ZW1iZXIsICIxMiI9RGVjZW1iZXIpICU+JQogIHNlbGVjdChEZWx0YSwgIjEiICwgIjIiICwgIjMiLCAiNCIsIjUiLCAiNiIsICI3IiwgIjgiLCAiOSIsICIxMCIsICIxMSIsICIxMiIpICU+JQogIHBpdm90X2xvbmdlcigtRGVsdGEsIG5hbWVzX3RvID0gIm1vbnRoIiwgdmFsdWVzX3RvID0gImRpc2NoYXJnZSIpCgojZmluZCBtYXggR1JEQyBtb250aCBmb3IgZWFjaCBkZWx0YQpEZWx0YU1heERpc2NoYXJnZSA8LSAKICBEZWx0YXNEaXNjaGFyZ2UgJT4lIAogIGdyb3VwX2J5KERlbHRhKSAlPiUgCiAgc2xpY2Uod2hpY2gubWF4KGRpc2NoYXJnZSkpICU+JSAKICByZW5hbWUoTWF4RGlzY2hhcmdlTW9udGggPSBtb250aCwgTWF4RGlzY2hhcmdlID0gZGlzY2hhcmdlKSAKCgpEZWx0YU1heERpc2NoYXJnZSRNYXhEaXNjaGFyZ2VNb250aCA9IGFzLm51bWVyaWMoRGVsdGFNYXhEaXNjaGFyZ2UkTWF4RGlzY2hhcmdlTW9udGgpCgoKI2pvaW4gd2l0aCBvdGhlciBkZWx0YSBkYXRhCkRlbHRhTWF4TWluRGlzY2hhcmdlIDwtIGxlZnRfam9pbihEZWx0YU1heE1pbiwgRGVsdGFNYXhEaXNjaGFyZ2UsIGJ5ID0gJ0RlbHRhJykKCiNjYWxjdWxhdGUgb2Zmc2V0CkRlbHRhTWF4TWluRGlzY2hhcmdlIDwtIERlbHRhTWF4TWluRGlzY2hhcmdlICU+JQogIG11dGF0ZSggRGlzc09mZiA9IGlmX2Vsc2UoYWJzKE1heERpc2NoYXJnZU1vbnRoIC0gTWF4TWVhbk5EU1NJbW9udGgpID4gNiwKICAgICAgICAgICAgICAgICAgICAgICAgICAgIDEyIC0gYWJzKE1heERpc2NoYXJnZU1vbnRoIC0gTWF4TWVhbk5EU1NJbW9udGgpLAogICAgICAgICAgICAgICAgICAgICAgICAgICAgYWJzKE1heERpc2NoYXJnZU1vbnRoIC0gTWF4TWVhbk5EU1NJbW9udGgpKQogICAgICAgICAgKQoKI0NvbXBhcmUgb2Zmc2V0IHdpdGggTkRTU0kgKGRlbHRhbWF4bWluJG1heG1lYW5ORFNTSW1vbnRoKQpnZ3Bsb3QoRGVsdGFNYXhNaW5EaXNjaGFyZ2UsIGFlcyh5ID0gRGVsdGEsIHggPSBEaXNzT2ZmKSkgKyBnZW9tX3BvaW50KCkgKyAKICBzY2FsZV94X2Rpc2NyZXRlKGxpbWl0cyA9IGMoMTo2KSwgYnJlYWtzID0gYygxOjYpKSArCiAgZXhwYW5kX2xpbWl0cyh4ID0gYygwLDYpKSAgKyAKICBnZ3RpdGxlKCJEaXNPZmZzZXQiKQpgYGAKCkxvb2sgYXQgR1JEQyBkYXRhIGJ5IGxhdGl0dWRlOgoKYGBge3J9CiNqb2luIGxhdCBkYXRhCkRlbHRhRGF0YXdMb2NhdGlvbnMgPC0gbGVmdF9qb2luKERlbHRhTWF4TWluRGlzY2hhcmdlLCBEZWx0YUxvY2F0aW9ucywgYnkgPSBjKCJEZWx0YSIgPSAiRGVsdGFzIikpCgpEZWx0YURhdGF3TG9jYXRpb25zIDwtIERlbHRhRGF0YXdMb2NhdGlvbnMgJT4lCiAgbXV0YXRlKEFic29sdXRlX0xhdGl0dWRlPSBhYnMoTGF0KSkKCiMgcGxvdCBvZmZzZXQgb24gZ3JhcGggYnkgbGF0CmdncGxvdChEZWx0YURhdGF3TG9jYXRpb25zLCBhZXMoeCA9IEFic29sdXRlX0xhdGl0dWRlLCB5ID0gRGlzc09mZikpICsgZ2VvbV9wb2ludCgpICsKICBzY2FsZV9jb2xvcl9ncmFkaWVudChsb3cgPSAieWVsbG93IiwgaGlnaCA9ICJyZWQiLCBuYS52YWx1ZSA9IE5BKSAKICAjKyBnZW9tX3Ntb290aChtYXBwaW5nID0gYWVzKHggPSBBYnNvbHV0ZV9MYXRpdHVkZSwgeSA9IERpc3NPZmYsICksIG1ldGhvZD1sbSApIAoKI2dnc2F2ZSgiRGlzT2Zmc2V0LnBkZiIsIHdpZHRoID0gNiwgaGVpZ2h0ID0gNCkKYGBgCgpBbmQgbm93IG9uIGEgbWFwOgpgYGB7cn0KI3Bsb3Qgb2Zmc2V0IG9uIG1hcApEZWx0YURpc09mZnNldE1hcCA8LSB3b3JsZCArCiAgZ2VvbV9wb2ludChhZXMoeCA9IExvbiwgeSA9IExhdCwgY29sb3IgPSBEaXNzT2ZmKSwKICAgICAgICAgICAgIGRhdGEgPSBEZWx0YURhdGF3TG9jYXRpb25zLCAKICAgICAgICAgICAgIHNpemUgPSA1KSArIHNjYWxlX2NvbG9yX2dyYWRpZW50KCBoaWdoID0gInJlZCIsIGxvdyAgPSAieWVsbG93IikgKwogIGdndGl0bGUoIk9mZnNldCBCZXR3ZWVuIEdSREMgZGlzY2hhcmdlIHBlYWsgYW5kIE5EU1NJIHBlYWsgaW4gd2F0ZXIiKQoKRGVsdGFEaXNPZmZzZXRNYXAKI2dnc2F2ZSgiRGVsdGFPZmZzZXRNYXAucGRmIiwgd2lkdGggPSA2LCBoZWlnaHQgPSA0KQoKYGBgCgoKCg==